home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 41 / Amiga Format CD41 (1999-06)(Future Publishing)(GB)[!][issue 1999-07].iso / -seriously_amiga- / programming / other / scm / slib / modular.scm < prev    next >
Text File  |  1999-04-19  |  5KB  |  159 lines
  1. ;;;; "modular.scm", modular fixnum arithmetic for Scheme
  2. ;;; Copyright (C) 1991, 1993, 1995 Aubrey Jaffer.
  3. ;
  4. ;Permission to copy this software, to redistribute it, and to use it
  5. ;for any purpose is granted, subject to the following restrictions and
  6. ;understandings.
  7. ;
  8. ;1.  Any copy made of this software must include this copyright notice
  9. ;in full.
  10. ;
  11. ;2.  I have made no warrantee or representation that the operation of
  12. ;this software will be error-free, and I am under no obligation to
  13. ;provide any services, by way of maintenance, update, or otherwise.
  14. ;
  15. ;3.  In conjunction with products arising from the use of this
  16. ;material, there shall be no use of my name in any advertising,
  17. ;promotional, or sales literature without prior written consent in
  18. ;each case.
  19.  
  20. (define (symmetric:modulus n)
  21.   (cond ((or (not (number? n)) (not (positive? n)) (even? n))
  22.      (slib:error 'symmetric:modulus n))
  23.     (else (quotient (+ -1 n) -2))))
  24.  
  25. (define (modulus->integer m)
  26.   (cond ((negative? m) (- 1 m m))
  27.     ((zero? m) #f)
  28.     (else m)))
  29.  
  30. (define (modular:normalize m k)
  31.   (cond ((positive? m) (modulo k m))
  32.     ((zero? m) k)
  33.     ((<= m k (- m)) k)
  34.     ((or (provided? 'bignum)
  35.          (<= m (quotient (+ -1 most-positive-fixnum) 2)))
  36.      (let* ((pm (+ 1 (* -2 m)))
  37.         (s (modulo k pm)))
  38.        (if (<= s (- m)) s (- s pm))))
  39.     ((positive? k) (+ (+ (+ k -1) m) m))
  40.     (else  (- (- (+ k 1) m) m))))
  41.  
  42. ;;;; NOTE: The rest of these functions assume normalized arguments!
  43.  
  44. (require 'logical)
  45.  
  46. (define (modular:extended-euclid x y)
  47.   (define q 0)
  48.   (do ((r0 x r1) (r1 y (remainder r0 r1))
  49.        (u0 1 u1) (u1 0 (- u0 (* q u1)))
  50.        (v0 0 v1) (v1 1 (- v0 (* q v1))))
  51.       ;; (assert (= r0 (+ (* u0 x) (* v0 y))))
  52.       ;; (assert (= r1 (+ (* u1 x) (* v1 y))))
  53.       ((zero? r1) (list r0 u0 v0))
  54.     (set! q (quotient r0 r1))))
  55.  
  56. (define (modular:invertable? m a)
  57.   (eqv? 1 (gcd (or (modulus->integer m) 0) a)))
  58.  
  59. (define (modular:invert m a)
  60.   (cond ((eqv? 1 (abs a)) a)        ; unit
  61.     (else
  62.      (let ((pm (modulus->integer m)))
  63.        (cond
  64.         (pm
  65.          (let ((d (modular:extended-euclid (modular:normalize pm a) pm)))
  66.            (if (= 1 (car d))
  67.            (modular:normalize m (cadr d))
  68.            (slib:error 'modular:invert "can't invert" m a))))
  69.         (else (slib:error 'modular:invert "can't invert" m a)))))))
  70.  
  71. (define (modular:negate m a)
  72.   (if (zero? a) 0
  73.       (if (negative? m) (- a)
  74.       (- m a))))
  75.  
  76. ;;; Being careful about overflow here
  77. (define (modular:+ m a b)
  78.   (cond ((positive? m)
  79.      (modulo (+ (- a m) b) m))
  80.     ((zero? m) (+ a b))
  81.     ((negative? a)
  82.      (if (negative? b)
  83.          (let ((s (+ (- a m) b)))
  84.            (if (negative? s)
  85.            (- s -1 m)
  86.            (+ s m)))
  87.          (+ a b)))
  88.     ((negative? b) (+ a b))
  89.     (else (let ((s (+ (+ a m) b)))
  90.         (if (positive? s)
  91.             (+ s -1 m)
  92.             (- s m))))))
  93.  
  94. (define (modular:- m a b)
  95.   (cond ((positive? m) (modulo (- a b) m))
  96.     ((zero? m) (- a b))
  97.     (else (modular:+ m a (- b)))))
  98.  
  99. ;;; See: L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  100. ;;; with Splitting Facilities." ACM Transactions on Mathematical
  101. ;;; Software, 17:98-111 (1991)
  102.  
  103. ;;; modular:r = 2**((nb-2)/2) where nb = number of bits in a word.
  104. (define modular:r
  105.   (ash 1 (quotient (integer-length most-positive-fixnum) 2)))
  106. (define modular:*
  107.   (if (provided? 'bignum)
  108.       (lambda (m a b)
  109.     (cond ((zero? m) (* a b))
  110.           ((positive? m) (modulo (* a b) m))
  111.           (else (modular:normalize m (* a b)))))
  112.       (lambda (m a b)
  113.     (let ((a0 a)
  114.           (p 0))
  115.       (cond
  116.        ((zero? m) (* a b))
  117.        ((negative? m)
  118.         "This doesn't work for the full range of modulus M;"
  119.         "Someone please create or convert the following"
  120.         "algorighm to work with symmetric representation"
  121.         (modular:normalize m (* a b)))
  122.        (else
  123.         (cond
  124.          ((< a modular:r))
  125.          ((< b modular:r) (set! a b) (set! b a0) (set! a0 a))
  126.          (else
  127.           (set! a0 (modulo a modular:r))
  128.           (let ((a1 (quotient a modular:r))
  129.             (qh (quotient m modular:r))
  130.             (rh (modulo m modular:r)))
  131.         (cond ((>= a1 modular:r)
  132.                (set! a1 (- a1 modular:r))
  133.                (set! p (modulo (- (* modular:r (modulo b qh))
  134.                       (* (quotient b qh) rh)) m))))
  135.         (cond ((not (zero? a1))
  136.                (let ((q (quotient m a1)))
  137.              (set! p (- p (* (quotient b q) (modulo m a1))))
  138.              (set! p (modulo (+ (if (positive? p) (- p m) p)
  139.                         (* a1 (modulo b q))) m)))))
  140.         (set! p (modulo (- (* modular:r (modulo p qh))
  141.                    (* (quotient p qh) rh)) m)))))
  142.         (if (zero? a0)
  143.         p
  144.         (let ((q (quotient m a0)))
  145.           (set! p (- p (* (quotient b q) (modulo m a0))))
  146.           (modulo (+ (if (positive? p) (- p m) p)
  147.                  (* a0 (modulo b q))) m)))))))))
  148.  
  149. (define (modular:expt m a b)
  150.   (cond ((= a 1) 1)
  151.     ((= a (- m 1)) (if (odd? b) a 1))
  152.     ((zero? a) 0)
  153.     ((zero? m) (integer-expt a b))
  154.     (else
  155.      (logical:ipow-by-squaring a b 1
  156.                    (lambda (c d) (modular:* m c d))))))
  157.  
  158. (define extended-euclid modular:extended-euclid)
  159.